Nonequilibrium transport through a spinful quantum dot with superconducting leads 
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We study the nonlinear cotunneling current through a spinful quantum dot contacted by two 
superconducting leads. Applying a general nonequilibrium Green function formalism to an effective 
Kondo model, we study the rich variation in the IV-characteristics with varying asymmetry in the 
tunnel coupling to source and drain electrodes. The current is found to be carried respectively by 
multiple Andreev reflections in the symmetric limit, and by spin-induced Yu-Shiba-Russinov bound 
states in the strongly asymmetric limit. The interplay between these two mechanisms leads to 
qualitatively different IV-characteristics in the cross-over regime of intermediate symmetry, consis- 
tent with recent experimental observations of negative differential conductance and re-positioned 
conductance peaks in sub-gap cotunneling spectroscopy. 
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During the last decade, a number of experiments 
have explored the problem of coherent electron trans- 
port through tunnel junctions and quantum dots cou- 
pled to superconducting leads [THZ] and the observation 
of sub-gap structure with steps near eV = 2A/n (n = 
1,2,3,...) in the low-bias current- volt age (IV) charac- 
teristics has been rationalized in terms of multiple An- 
dreev reflections (MAR) [8- 10 . Interestingly, recent ex- 
periments have revealed that odd-occupied (Coulomb 
blockaded) quantum dots exhibit novel transport char- 
acteristics not included within standard MAR- mediated 
transport [31 [TIHT5] . 

The physics of odd-occupied quantum dots coupled to 
superconducting leads is intimately tied to the study of 
magnetic impurities in superconductors. Magnetic impu- 
rities are known to generate localized Yu-Shiba-Rusinov 
(YSR) bound states inside the superconducting gap, off- 
set from the gap edge roughly by the magnitude of the 
exchange coupling[16 . These bound states serve as in- 
formative signatures in tunneling spectroscopy of super- 
conductor surfaces [T6l [T7] and very recently they have 
also been probed by scanning tunneling spectroscopy of 
MnPc molecules adsorbed on Pb(lll) surfaces [18 j. The 
question remains, however, which signatures are to be ex- 
pected in a generic voltage-biased quantum dot transport 
experiment, where the system is driven out of equilibrium 
by a finite Josephson frequency 2eV/h. 

Recent theoretical work studied Andreev bound states 
induced by a spinful quantum dot coupled to a supercon- 
ductor, with focus on the competition between Kondo 
and superconducting singlet correlations [19 -32J. From 
detailed calculations of low-energy spectral properties, 
one may thus infer about e.g. the Josephson current car- 
ried by bound states. Far less is known about nonlinear 
transport through an interacting, voltage-biased quan- 
tum dot, the focus of the present paper [33U36]. 

In this Letter, we perform a theoretical study of 
spin-induced YSR bound states and their influence on 



the nonequilibrium transport through S/QD/S cotunnel 
junctions. This study is motivated by recent experiments 
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FIG. 1: (Color online) The top figure shows the conductance 
versus gate (Vg a ^ e ) and bias (Vg^) voltage for a carbon nan- 
otube quantum dot coupled to superconducting leads reveal- 
ing a qualitatively different sub-gap structure at fillings 3, 
7, and 11. For more details on the experimental setup see 
Ref. [14] . The bottom figure shows calculated IV character- 
istics through a S/QD/S cotunnel junction obtained within 
the present approach. For all curves Jll/A = 11.9 and (top 
to bottom) Jrr/A = 0.05,0.8,3.2,11.9. The curves are off- 
set for clarity. Top inset: Sketch of the S/QD/S cotunnel 
junction consisting of a spinful quantum dot tunnel-coupled 
to the superconducting leads. Bottom inset: Second order 
self-energy diagram included in the calculation. 
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on quantum dots where the gate voltage can alter the 
nature of the quantum "impurity" as shown in Figjl] at 
specific odd fillings (3, 7, 11), the conductance through 
this carbon nanotube quantum dot is strikingly different 
from all the even-occupied fillings. The unusual behavior 
includes both negative differential conductance (NDC) 
and re-positioned conductance peaks inside the bias volt- 
age range of ±2A/e. Our modeling consists of a gen- 
eral framework for the nonequilibrium current through 
an interacting dot which is treated to second order in 
the lead-dot tunneling. Our treatment, therefore, does 
not include Kondo correlations which are, as we show 
below, not responsible for the unusual conductance fea- 
tures in the experimental plot in Figjl] Rather, as shown 
in the lower panel in Figjl] the low-bias Andreev trans- 
port is significantly altered by the presence of magnetic 
bound states which, in the asymmetric case, completely 
dominate the transport and generate NDC and shifted 
sub-gap conductance peaks similar to the experimental 
data. Specifically, the calculated IV curves in Figjl] are 
obtained for fixed (varying) coupling to the left (right) 
lead. The bottommost curve is for a symmetric junc- 
tion whereas the top curve is in the tunneling limit with 
highly asymmetric coupling. In the parameter regime 
between these two limiting cases, one clearly observes a 
transition from "conventional" MAR transport to a new 
double- step IV curve with associated NDC in the asym- 
metric limit. In addition, a regime of intermediate cou- 
pling asymmetry exists where mainly the current near 2 A 
is altered (see e.g. the red curve) reminiscent of fillings 
7 and 11 in the experimental conductance map of Figjl] 
The absence of NDC at odd fillings 5 and 9 can be ex- 
plained by a more symmetric coupling of the lower-lying 
nanotube orbital to the leads [14 . 

We consider a system of two superconducting leads 
bridged by a quantum dot with its highest partially oc- 
cupied orbital represented by a single-orbital Anderson 
model. By focussing on the cotunneling regime well in- 
side a Coulomb diamond, where charge fluctuations hap- 
pen only virtually, a Schrieffer- Wolff transformation [37] 
results in the following effective Hamiltonian 



ak 

+ I ^2 Joc' a S l ^ a , k ,m l xl) ak . 



(1) 



a' k' ,ak ,i 



Here, ^ fe = (c* fct , c a _ fct , c^, c a - k ^ , £ ak is the quasi- 

particle energy, A a = A' a + zA" is the BCS-gap in 
lead a = L,R, and the Nambu matrices are defined by 



m 



r 1 (g)r 2 , m 2 



)r 2 , m 3 



T 2 and m z 



= t°(8>t , m — r^r , 
t 3 , where r % denote the 



Pauli matrices and r° is the 2x2 identity matrix. S % is 
the quantum spin on the dot and J a 'a is the exchange- 
cotunneling amplitude with lead index a, ol = L,R, 



which satisfies the relation JlrJrl = JllJrr- 

The applied bias voltage, included through Hy = 
~ J2 a Ma^a with fjL a = ±V/2, may be incorporated 
into the spinors = ^L^K^™^ (»=!), 

where the time-dependence of ^ ak ^(t) is now deter- 
mined by the unbiased (non-interacting) leads. The cor- 
responding contour ordered Green function transforms as 
G ari , aW (t,t') = e i ^< i G a ^ v (M0e~ itV "' m -''>' ) and 
G can now be determined from perturbation theory in 
the transformed cotunneling term: 



ot'k'r)' ,akr],i 



with the time- dependent amplitudes J^'ait) = 
e"^' m ^'J^ Q e^<. The problem of voltage-biased 
superconducting junctions is inherently time-dependent 
and the Green functions acquire an explicit dependence 
on two times, t and t r . However, because of the harmonic 
dependence of the perturbation, the Green functions 
qR,a,>,< can k e ex p reS g ec i as a Fourier series in T = t+t f 



G(t,t') = G(T,t) =J2e inTV/2 Gn(r), 



(3) 



where r = t — t' is the relative time and G n (uj) = 
dre ZUJT G n (r), with G being shorthand for any of the 
functions G <,>,R,A . 

In the following, the Green functions and the self- 
energy are written as 8 x 8 matrices containing a = L,R 
and the four Nambu indices n = 1,...,4. Notice that 
the redundancy in the 4-spinor notation gives rise to 
anomalous Wick-contractions. Being identical replicas, 
however, these are effectively included by multiplying all 
self-energy diagrams by a factor of 4 and then^employing 
the normal Dyson equation written here for G r (uj)\ 



G r (uj) = 5 n0 G R W(uj) + G R(0 \co - nV/2) 
x X R -m(" ~ ™V/2)G r (uj + (n - m)V/2). 



(4) 



The solution generally takes the form 

G r (lo - nV/2) = [(1 - M R (uj))~ 1 \ nQ G R{ - \oj), (5) 

where M R m (uj) = G R( -°\u - nV)E R _ m (u; - (m + n)V/2) 
and the free (fc-summed) Green function is given by 



G r ^\uj) 



ujrn 



A m 1 + A m 2 



(6) 



using a constant density of states, vp = 1/(212), with 
D denoting the bandwidth. Using the irreducible self- 
energy, calculated to leading order in perturbation the- 
ory, Eq.([5| can be solved by numerical matrix inversion 
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FIG. 2: (Color online) IV curves for a S/QD/S cotunnel 
junction with symmetric couplings, Jll — Jrr — Jlr — 
J. From top to bottom the curves correspond to J/ A = 
22.6, 17.3, 12.7,8.8, 5.7. All curves are normalized by Gn de- 
termined for J/A — 22.6. Inset: dot spectral function for 
J I A = 12.7 with V = (solid) and V = A/e (dashed). 



FIG. 3: (Color online) IV curves for a S/QD/S co- 
tunnel junction with coupling asymmetry Jrr/ Jll = 
0.004,0.007,0.016,0.062,1.0 for lines 1-5, respectively. All 
curves have Gn — 0.013(2e 2 //i), corresponding to vfJll ~ 
0.744, 0.511, 0.335, 0.170, 0.042. Inset: as in Fig[2]for the case 
with Jrr/ Jll = 0.016. 



for each frequency uj after suitable truncation of the num- 
ber of included harmonics n. Having done this matrix 
inversion, one obtains the T-matrix components 



T^H = ^S«_ m ( w -m^/2)[(l-M( w + n^/2))- 1 ] 



m0 5 



(7) 



m 

+ G^<(oj)E A (Lo) m ]] [{l-M^ + nV^))- 1 ]^ (8) 
which enter the time-averaged current as 



(i) 



duj 



(9) 



This expression follows from the underlying Anderson 
model, and due to the_ energy dependence of the BCS 
density of states both T R and T < must be computed in 
order to determine the current. Given an expression for 
the conduction electron self-energy, this current formula 
applies to any dot Hamiltonian. The numerical cost lies 
in the matrix inversion, which becomes increasingly de- 
manding as bias voltage is reduced, and a larger number 
of harmonics must be included to ensure convergence. 

For even-occupied (spinless) dots, the first order self- 
energy is exact and formulas (|6][9| quantitatively repro- 
duce the IV curves presented, e.g. in Refs. [8HT0]. Hence- 
forth, we focus our attention to odd occupied dots where 
the spin turns out to play an important role. We restrict 
the discussion to the weak-coupling regime, vpj < 1 



(J = max{ Ja'a})-, and retain only the second order self- 
energy (see diagram in the lower inset of FigJI]): 

S 7) n( w ) = S J^a-Ja-arn^m^ (10) 

?7 / 77,i,Q; // 

x G^?\u - [ti a >m 3 YY - ^"ml,^ + nV/2]), 

together with the corresponding advanced and lesser 
components. For all numerical evaluations presented in 
Figs. [T][3| we take D = 20 A. 

For equal coupling to the two leads, the IV curves for 
the spinful dot are shown in Figj2j Evidently, the current 
is dominated by MAR steps at 2A/n, resembling that 
of an even occupied dot, and there appears to be little 
effect of the spin. This is in stark contrast to the case 
where tunnel couplings to the two leads are very different. 
As seen in Fig(3J the coupling asymmetry significantly 
enhances the sub-gap current and the IV characteristics 
now exhibit a marked cusped two-step shape. We can 
understand this evolution from symmetric to asymmetric 
junctions by taking a closer look at the physical retarded 
T-matrix, obtained as: 



where uj' = uj 



Figs. 



(/i a m^ + /i a /m^,)/2. In the insets of 



we plot — ImT£^ l1 .${uj + V/2)/tt, which is pro- 
portionaTto the time- averaged dot-electron spectral func- 
tion for the underlying Anderson model, for respectively 
V = and V = A/e. In the symmetric limit (Fig{2| 
sub-gap features are predominantly due to MAR, and for 
these disperse strongly with time and are therefore 
smeared out by time averaging. This results in broad res- 
onances which depend significantly on the bias voltage. 
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In the asymmetric limit (Fig j3J) the spectrum shows a pair 
of YSR Andreev bound states [16l [38] located symmetri- 
cally around the left chemical potential at hl±^s where 
ujs = A(l—x)/(l + x) with x = 3(ttuf Jll/^) 2 and hence 
crossing; each other when 

v f Jll = 4/(ttV3) « 0.74, cor- 
responding to curve 1 in Fig. [3] Note that the peaks 
in the inset have been smeared with an artificial broad- 
ening chosen small enough that it does not affect the 
calculated current. In contrast to the symmetric case, 
these bound states do not disperse with time and hence 
survive the time- averaging. This is seen explicitly in the 
inset of Fig. [3] where peaks at uos = 0.66A simply shift 
by eV/2 upon an applied bias voltage. The values of 
us are directly related to the cusps in the current at 
eV = A ± ojs(= 1.66,0.34) (case 3 in the main panel of 
Figj3|, corresponding to the matching up of YSR bound 
states with the coherence peaks of the weakly probing 
right lead. In this case, the current is dominated by An- 
dreev processes matching up with the symmetrically posi- 
tioned bound states. This is similar to the resonant MAR 
discussed for noninteracting junctions [39, 40 but now 
the resonant effect is enhanced because the bound states 
are positioned symmetrically around and therefore 
match both an incoming electron and an outgoing hole 
of a contributing Andreev reflection process. Finally, a 
strong coupling asymmetry is seen to change the quasi- 
particle tunneling step at V = 2A/e into a local current 
minimum, consistent with the experimental observations 
of NDC in Refs. [HI HI HI EH], and seen also in Figjl] 

The progression of curves in Fig |3] corresponds to in- 
creasing asymmetry factor Jll/Jrr with fixed normal 
state conductance, Gn — (37r 2 /4)|^ J^ j R| 2 (2e 2 /h) (cal- 
culated to second order), and therefore also to increas- 
ing values of Jll- For sufficiently weak couplings, the 
dominant higher order self-energy corrections may be in- 
cluded by using the renormalized exchange coupling ob- 
tained from Poor man's scaling from D down to A [41]. 
Thus replacing VfJll with l/ln(A/Tx), we conclude 
that the two YSR bound states lie at zero energy when 
l/ln(A/7V) « 0.74, i.e. T K « 0.26 A. This is consistent 
with numerical renormalization group calculations for the 
equilibrium problem[20, 23, 29 , which show a transition 
from doublet to singlet groundstate, corresponding to a 
crossing of the bound states, at Tk ~ 0.3A. For even 
stronger couplings, the YSR bound states again lie sym- 
metrically around hl and our calculations would give IV 
curves which look qualitatively similar to those in Figj3] 
For such couplings, however, the Kondo screening is so 
strong that we can no longer rely on our perturbative ap- 
proach. Note, however, that our approach would be exact 
in all regimes if dealing with a classical magnetic impu- 
rity, potentially relevant for experiments probing impu- 
rities on surfaces [TT] [18] . 

In summary, we have presented a general Hamilto- 
nian approach to the nonequilibrium current through an 
interacting quantum dot contacted by superconducting 



leads. For odd-occupied dots, the spin leads to IV charac- 
teristics being strongly dependent on coupling asymme- 
try, giving rise to negative differential conductance and 
shifted sub-gap conductance peaks in agreement with ex- 
periments. More work is required to study the asymme- 
try variations in the cross-over regime towards stronger 
couplings, Tk > A, where Kondo screening eventually 
quenches local superconducting fluctuations. 
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